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ABSTRACT 

This  report  describes  the  development  of  a  two-dimensional  multi-material  Eulerian 
hydrocode  to  model  the  effects  of  detonating  condensed  phase  explosives  on 
surrounding  materials.  The  code  solves  the  Euler  equations  for  the  conservation  of 
mass,  momentum,  and  energy  for  an  inviscid,  compressible  fluid.  Operator  splitting  is 
used  to  reduce  the  two-dimensional  calculation  into  coupled  one-dimensional 
equations,  which  are  then  solved  using  the  Flux-Corrected  Transport  (FCT)  algorithm 
of  Boris  and  Book.  Non-reacting  materials  are  described  using  either  a  perfect  gas,  Mie- 
Gruneisen,  or  Tait  equation  of  state,  while  the  energetic  materials  are  described  using 
either  a  BKW  equation  of  state  and  Forest  Fire  reaction  rate  model,  or  the  JWL  equation 
of  state  and  the  Ignition  and  Growth  reaction  rate  model.  A  modified  Young's 
algorithm  is  used  to  maintain  a  sharp  interface  between  different  materials  on  die 
computational  mesh.  A  brief  description  of  the  major  components  of  the  coding  is 
provided,  and  then  several  applications  of  the  code  are  described,  including  the 
simulation  of  bullet  impact  experiments,  the  underwater  sympathetic  detonation  test, 
and  the  modified  gap  test. 
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Numerical  Simulation  of  Detonation  in 
Condensed  Phase  Explosives 


Executive  Summary 

The  ability  to  numerically  simulate  the  propagation  or  failure  of  detonation  in  a  variety 
of  gaseous  and  condensed  phase  explosives  is  of  fundamental  interest  to  the  military 
community.  Much  of  the  work  of  the  DSTO  is  now  done  through  computer  simulation 
techniques,  and  hydrocodes  are  die  favoured  method  for  modelling  warhead  outputs 
and  their  effects  on  targets.  Enhancing  the  capabilities  of  our  hydrocodes  directly 
enhances  our  ability  to  provide  better  advice  to  the  Australian  Defence  Force  in  die 
areas  of  airblast,  missile  defeat,  fragmentation,  and  the  defeat  of  fixed  or  mobile  land 
targets.  Hydrocodes  containing  advanced  energetic  reaction  rate  models  are  also  used 
to  model  safety  related  aspects  of  weapons  such  as  sympathetic  detonation,  cookoff, 
and  bullet  and  fragment  impact.  The  work  described  here  has  provided  WSD  with  a 
sophisticated  hydrocode  which  can  be  used  to  provide  an  in-house  capability  to  model 
non-standard  explosive/ shock  interaction  problems. 

In  this  report  we  describe  the  development  of  an  in-house  two-dimensional  multi¬ 
material  Eulerian  hydrocode  to  model  the  effects  of  detonating  condensed  phase 
explosives  on  surrounding  materials.  The  code  has  the  capability  to  model  standard 
military  explosives,  which  have  so  called  ideal  detonation  behaviour,  as  well  as  the 
newer  underwater  military  explosives  which  display  highly  non-ideal  behaviour.  The 
code  has  been  designed  to  provide  DSTO  personnel  involved  in  the  modelling  of 
energetic  materials  response  with  an  in-house  hydrocode  capability  which  can  be 
quickly  adapted  to  model  novel  materials  in  non-standard  environments  for  which 
current  commercial  hydrocodes  may  not  be  appropriate. 

The  code  solves  a  set  of  equations  describing  the  conservation  of  mass,  momentum, 
and  energy  for  an  inviscid,  compressible  fluid.  A  computational  technique  known  as 
operator  splitting  is  used  to  reduce  the  two-dimensional  calculation  into  sets  of 
coupled  one-dimensional  equations,  which  are  then  solved  using  the  Flux-Corrected 
Transport  (FCT)  algorithm  of  Boris  and  Book.  Non-reacting  materials  are  described 
using  either  a  perfect  gas,  Mie-Gruneisen,  or  Tait  equation  of  state.  Ideal  explosives  are 
described  using  either  a  BKW  equation  of  state  and  Forest  Fire  reaction  rate  model,  or 
the  JWL  equation  of  state  and  the  Ignition  and  Growth  reaction  rate  model.  A 
programmed  bum  option  is  also  available,  which  forces  the  explosive  to  detonate  at 
the  correct  detonation  velocity  and  pressure.  A  model  for  non-ideal  explosives  is  also 
included.  This  uses  a  three  term  reaction  rate  model  developed  for  composite 
explosives,  and  a  polytropic  equation  of  state  with  a  density  dependent  index. 

The  code  described  here  is  an  Eulerian  code,  which  means  that  the  computational  grid 
is  fixed  in  space  and  different  materials  in  the  simulation  flow  through  the  grid  at 
different  rates.  Maintaining  a  sharp  interface  between  different  materials  is  a  common 
problem  in  Eulerian  codes,  and  the  problem  is  solved  here  by  using  a  modified 
Young's  algorithm  to  prevent  numerical  diffusion  from  smearing  the  interface  between 
different  materials  on  the  computational  mesh.  This  is  a  two-dimensional  scheme  in 
which  an  interface  between  two  different  materials  within  a  cell  is  approximated  by  a 
straight  line. 
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1.  Introduction 


The  ability  to  numerically  simulate  the  propagation  or  failure  of  detonation  in  a  variety 
of  gaseous  and  condensed  phase  explosives  is  of  fundamental  interest  to  the  military 
community.  Computer  codes  which  allow  these  calculations  to  be  performed  in  a 
routine  manner  were  developed  during  the  1960’s  and  1970's,  and  many  of  these  codes 
have  now  been  released  to  the  general  defence  community.  The  codes  developed  by 
Mader  at  the  Los  Alamos  National  Laboratory  to  model  detonation  in  condensed  phase 
explosives  include  the  Lagrangian  finite  difference  codes  SIN  [1]  and  2DL  [2],  as  well 
as  the  multimaterial  Eulerian  codes  2DE  [3]  and  3DE  [4].  An  extensive  review  of  the 
theory  and  equations  used  in  these  codes  has  been  given  in  the  monograph  by  Mader 
[5].  Other  codes  have  been  developed  which  concentrate  primarily  on  material 
deformation  and  unreactive  shocks,  and  these  include  the  HELP  code  of  Hageman  and 
Walsh  [6],  the  HULL  code  [7],  and  the  DYNA2D  code  of  Hallquist  [8].  These  codes 
contain  only  a  rudimentary  treatment  of  the  physics  of  detonation  when  compared 
with  Mader's  codes,  although  recent  versions  of  DYNA2D  have  included  the  Lee  and 
Tarver  Ignition  and  Growth  model  for  the  development  of  detonation  in  condensed 
phase  explosives . 

Over  the  past  12  years  Weapons  Systems  Division  has  acquired  most  of  these  codes 
and  applied  them  to  a  variety  of  explosive /material  interaction  problems.  These  have 
included  the  use  of  HELP  to  model  jet  formation  in  the  MRL  38  mm  shaped  charge  [9], 
the  use  of  SIN  and  2DL  to  model  the  propagation  of  detonation  in  ladder  fracture  tape 
[10],  and  the  application  of  DYNA2D  to  model  liner  deformation  in  a  variety  of 
explosively  formed  projectiles  [11-14].  While  these  examples  have  been  notable 
successes,  there  have  also  been  problems  when  codes  have  been  applied  to  situations 
for  which  they  were  initially  unsuited.  Lack  of  documentation  regarding  the  inner 
working  of  the  codes  has  also  often  caused  problems,  and  the  large  size  of  some  of 
these  multipurpose  codes  has  also  been  a  disadvantage  when  codes  have  had  to  be 
migrated  between  different  computing  systems. 

To  overcome  some  of  these  problems  WSD  has  been  developing  in-house  reactive 
Eulerian  hydrocodes  using  the  Flux-Corrected  Transport  (FCT)  algorithm  of  Boris  and 
Book  [15,16].  Early  work  centred  on  the  development  of  a  multi-species  code  to  model 
detonation  transfer  in  gaseous  explosives  [17-19],  and  the  development  of  codes  to  model 
shock  propagation  in  elastic-plastic  materials  [20,  21].  More  recent  applications  of  these 
codes  have  included  the  modelling  of  detonation  initiation  and  failure  in  gaseous 
systems  [22-26],  and  air  blast  problems  associated  with  explosive  breaching  operations 
[27-31].  The  codes  have  recently  been  considerably  extended  by  the  addition  of  advanced 
interface  tracking  algorithms  [32,33],  enhanced  graphics  capabilities  [34],  and  the 
inclusion  of  algorithms  to  model  detonation  in  condensed  phase  explosives  [35]. 

The  code  development  described  above  has  resulted  in  the  production  of  an  in-house 
2D  multi-material  reactive  Eulerian  hydrocode  (known  as  MULTI),  which  has  recently 
been  used  to  model  several  explosive /material  interaction  problems  concerned  with 
explosive  sensitivity  and  safety  [36,  37].  In  this  report  we  describe  the  basic  features  of 
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the  MULTI  code  and  give  examples  of  its  application  to  several  standard  detonation 
problems.  MULTI  has  the  capability  to  model  standard  military  explosives,  which  have 
so  called  ideal  detonation  behaviour,  as  well  as  the  newer  underwater  military 
explosives  which  display  highly  non-ideal  behaviour. 

In  the  next  section  we  describe  the  basic  hydrodynamic  transport  equations  which 
govern  the  propagation  of  shocks  in  both  reactive  and  non  reactive  materials.  In  section 
3  we  discuss  the  modelling  of  ideal  explosives  and  describe  their  behaviour  using  both 
the  Becker-Kistiakowsky-Wilson  (BKW)  equation  of  state  and  the  Forest  Fire  reaction 
rate  model,  as  well  as  the  Jones-Wilkins-Lee  (JWL)  equation  of  state  and  the  Ignition 
and  Growth  reaction  rate  model.  In  section  4  we  consider  non-ideal  explosives  and 
describe  a  model  for  the  reaction  rate  and  equation  of  state  of  a  typical  underwater 
explosive,  PBXN-111.  The  numerical  solution  of  the  coupled  set  of  equations  describing 
hydrodynamic  flow  and  chemical  reaction  is  described  in  section  5.  We  first  briefly 
discuss  the  development  of  MULTI,  and  then  give  a  broad  outline  of  the  function  of 
each  of  the  various  files  and  subroutines  which  are  contained  in  the  code,  as  well  as  an 
explanation  of  the  input  data  file.  In  section  6  we  describe  the  results  of  the  code  when 
used  to  simulate  typical  problems  for  both  ideal  and  non-ideal  explosives,  and  in 
section  7  we  provide  a  summary  of  the  code  capabilities  and  areas  of  application. 


2.  HYDRODYNAMIC  FLOW  EQUATIONS 

2.1  Transport  Equations 

Detonation  phenomena  occur  on  the  microsecond  time  scale,  and  it  is  usual  to  assume 
that  energy  transport  by  heat  conduction,  viscosity,  and  radiation  is  negligibly  small 
compared  with  transport  by  motion.  The  pressures  generated  by  the  detonation  of  a 
condensed  phase  explosive  are  of  the  order  of  105  atmospheres.  Under  these  conditions 
the  strength  of  any  confining  material  is  completely  negligible,  and  the  material 
responds  hydrodynamically.  In  this  case,  the  appropriate  equations  which  describe  the 
material  response  are  the  Euler  equations,  which  describe  the  conservation  of  mass, 
momentum  and  energy  for  an  inviscid,  compressible  fluid.  These  have  the  following 
form  [16]: 

^+v-Gw)  =  o  (i) 

Ot 

^Uv.(,cw)  =  ~VP  (2) 

Ot 

8F 

—  +  V.(£v)  =  -V.(Pv)  (3) 

ot 

In  equations  (1)  -  (3),  p  is  the  density,  v  the  fluid  velocity,  P  the  hydrodynamic 
pressure,  and  E  the  total  energy  per  unit  volume,  which  is  given  by  the  following 
expression 
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E  =  pe  +  0.5pv2  (4) 

where  e  is  the  specific  internal  energy. 

The  above  equations  are  supplemented  by  both  a  mechanical  equation  of  state  for  the 
pressure  and  a  caloric  equation  of  state  for  the  temperature,  which  have  the  form 

P  =  P  (p,  e.  A.)  (5) 

T  =  T  (p,  e,  X)  (6) 

where  T  is  the  temperature  and  X  is  the  explosive  mass  fraction.  The  equations  are  then 
closed  by  specifying  an  expression  for  the  rate  of  change  of  X.  The  exact  form  of  this 
expression  depends  on  the  details  of  the  explosive  decomposition  reaction,  and  for  the 
moment  we  simply  write 

^  =  /(P,e,2)  (7) 

at 

When  the  reaction  occurs  in  a  fluid  moving  with  a  velocity  v  the  time  derivative  in 
equation  (7)  denotes  a  derivative  following  the  fluid  flow,  and  equation  (7)  becomes 

^  +  v-VA  =  /(P,e„l)  (8) 

Ot 

Equation  (8)  can  then  be  combined  with  equation  (1)  so  that  the  equation  describing 
chemical  reaction  in  the  moving  fluid  can  be  written  in  the  form 

^l  +  V.(p2v)=/(F,e,A)  (9) 

Ot 

Equation  (9)  then  has  the  same  form  as  equations  (1)  through  (3)  and  can  be  solved  by 
the  same  numerical  technique. 

The  numerical  solution  of  the  progress  of  the  detonation  therefore  requires  the  solution 
of  equations  (1),  (2),  (3),  which  describe  conservation  of  mass,  momentum  and  energy 
respectively,  together  with  the  solution  of  the  equations  of  state  given  by  equations  (5) 
and  (6),  and  the  equation  describing  the  rate  of  reaction,  equation  (9).  The  numerical 
techniques  used  to  solve  these  equations  will  be  described  in  section  5.  In  the  next  few 
sections  we  discuss  appropriate  equations  of  state  which  are  used  to  describe  the 
behaviour  of  both  non-energetic  and  reactive  materials. 

2.2  Equations  of  state  for  non-reactive  materials 

The  most  commonly  used  equation  of  state  for  condensed  phase  materials  is  the  Mie- 
Gruneisen  equation  of  state,  which  has  the  form  [38]: 

P  =  PH  +  y(e-eH)p  (10) 

where  Ph  and  eH  refer  to  the  pressure  and  specific  internal  energy  along  the  shock 
compression  curve,  known  as  the  Hugoniot,  and  which  have  the  form 
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Ph  =  Po  +  poCV)/  (l-s0ri)2  (11) 

eH  =  e0  +  0.5(Ph  +P0)r|/  P  o  (12) 

where  r|  is  the  compression,  defined  by  r\  -  1  -  p0/p,  and  cD  is  the  sound  speed  in  the 
undisturbed  material.  In  the  above  expressions  the  shock  velocity  us  and  the  particle 
velocity  up  are  assumed  to  obey  a  linear  relationship  of  the  form  us  =  c0  +  s0up.  The 
constant  y  is  known  as  the  Gruneisen  gamma,  and  the  product  yp  is  assumed  to  be  a 
constant,  and  is  given  by  the  expression 

yp  =  3a  /  (CVK)  (13) 

where  a  is  the  coefficient  of  linear  thermal  expansion,  K  is  the  isothermal 
compressibility,  and  Cv  is  the  specific  heat  at  constant  volume. 

Equation  (10)  effectively  describes  the  hydrodynamic  shock  response  of  many 
materials,  including  metals,  unreacted  explosives,  and  water.  In  the  case  of  water 
however  a  simpler  equation  of  state  ,  the  Tait  equation  [39],  is  sometimes  used.  This 
has  the  form 


p_  = 

f p ' 

n 

-1 

+ 1  for  —  >1 

Po 

\Po) 

\  Po  J 

Po 

P_ 

Po 


KPo) 


for  -£-<a 
Po 


(15) 


where  po  and  po  are  pressure  and  density  at  standard  conditions,  taken  as  1  atm  and 
1000  kg/m3,  respectively.  The  quantities  B,  n,  and  m  are  constants  which  have  the 
values  3010  atm,  7.14,  and  21457.15  respectively  [40]. 


3.  IDEAL  EXPLOSIVES 

3.1  The  HOM  Equation  of  state 

Energetic  materials  require  different  equations  of  state  to  describe  the  condensed  phase 
material  before  detonation,  and  the  gaseous  detonation  products  after  reaction.  The 
HOM  Equation  of  State  (EOS)  uses  the  Mie-Gruneisen  EOS  to  describe  the  unreacted 
material,  while  the  BKW  EOS  is  used  to  describe  the  highly  non-ideal  gaseous  products 
after  reaction.  For  a  mixture  of  gases  the  BKW  EOS  has  the  form  [5]: 

PgVs/(RT)  =  l  +  x*exp([3x)  (16) 

where  x  =  k/[vg(T+9)“]  and  k  is  the  average  covolume  defined  in  terms  of  the 
individual  covolumes  ki  as  k  =  kZx;  k„  where  Xi  is  the  mole  fraction  of  species  I.  The  a. 
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p,  k,  and  0  are  constants  adjusted  to  reproduce  the  detonation  pressure  and  velocity 
obtained  experimentally. 

By  using  various  thermodynamic  expressions  equation  (16)  can  be  used  to  construct 
the  free  energy  for  a  given  composition,  and  then  the  equilibrium  composition  can  be 
determined  by  minimizing  the  free  energy  with  respect  to  the  molar  fractions.  This 
procedure  is  described  in  detail  in  the  monograph  by  Mader  [5],  and  in  a  recent  DSTO 
report  [41].  To  avoid  performing  this  calculation  every  time  the  equation  of  state  of  the 
products  is  needed  in  a  hydrocode  calculation  however  it  is  assumed  that  the  states 
encountered  in  a  typical  calculation  are  close  to  those  on  the  eqilibrium  isentrope 
through  the  CJ  point.  Mader  has  used  his  BKW  code  to  fit  the  equilibrium  isentrope 
through  the  CJ  point  for  many  different  explosives  to  a  set  of  polynomials  of  the  form 

In  JFf  =  A  +  B(lnvg)  +  C(lnvg)2  +  D(  lnvg)3  +  E(in  vg)4  (17) 

In*:  =K  +  L(kiPi)  +  M(lnPi)2+N(lnfy3+0(lnPi)4  (18) 

ln7]  =  <2  +  i?(lnvg)  +  S(lnvg)2  +  r(lnvg)3  +  £/(lnvg)4  (19) 

where  vs  is  the  specific  volume  of  the  detonation  products,  the  shifted  internal  energy 
s\  =  si  +  Z ,  and  Z  is  a  constant  used  to  change  the  standard  state  to  be  consistent  with 
the  condensed  explosive  one.  In  a  hydrocode  calculation  then  the  BKW  equation  of 
state  for  the  reaction  products  is  approximated  by  an  equation  of  Mie  Gruneisen  form, 
where  the  reference  state  is  now  the  isentrope  through  the  CJ  point,  rather  than  the 
Hugoniot,  ie. 

Ps  =  Pi  +  y  p  (eg  -  ei)  (20) 

where  the  subscripts  i  and  g  refer  to  the  isentrope  and  the  gaseous  products,  and  y  is 
now  the  Gruneisen  gamma  for  the  detonation  products. 

3.2  The  JWL  Equation  of  state 

The  Jones-Wilkins-Lee  equation  of  state  has  been  used  extensively  to  describe  the 
behaviour  of  the  detonation  products  of  explosives  in  contact  with  metal  casings.  It  has 
the  form  [42]: 

P  =  A 1  -  ■§=)  exp  H,  V)  +  B(l—E-)  exp  (-R2F)  +  ~  (21) 

iV|  r  Jt\2  r  r 

where  p  is  the  pressure  in  megabars.  V  and  E  are  given  by  V  =  p0/p  and  E  =  pce,  and  A, 
B,  Ri,  R2,  and  co  are  constants  for  a  particular  explosive.  The  JWL  constants  are  derived 
by  a  trial  and  error  procedure  in  which  the  expansion  record  in  a  standard  cylinder 
expansion  test  is  simulated  using  a  hydrocode,  usually  DYNA2D,  and  the  constants  are 
varied  until  the  simulated  expansion  record  matches  the  experimental  data.  Values  for 
these  constants  can  be  found  in  a  variety  of  sources,  including  Dobratz  and  Crawford 
[43],  and  Karpp  [44].  Equation  (21)  can  also  been  used  to  describe  the  equation  of  state 
of  unreacted  explosives  [51]. 
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3.3  Reaction  rate  models 


The  shock  initiation  of  a  condensed  phase  heterogeneous  explosive  is  a  complicated 
process  which  is  still  not  completely  understood.  Most  condensed  phase  explosives 
consist  of  poly  crystalline  materials  containing  voids  of  various  shapes  and  sizes,  defect 
structures,  and  often  small  amounts  of  polymeric  binders  and  plasticisers.  When  a 
shock  wave  travels  through  such  material  it  provides  heating  both  by  bulk 
compression  and  by  the  interaction  of  the  shock  with  the  various  density 
discontinuities  and  defect  structures.  The  localised  regions  of  high  temperature  caused 
by  these  density  discontinuities  are  known  as  hot  spots,  and  if  conditions  are 
favourable  these  hot  spots  may  begin  to  react  and  lead  to  the  formation  of  a  stable 
detonation,  even  though  the  temperature  rise  caused  by  the  bulk  heating  may  be 
insufficient  by  itself  to  lead  to  detonation.  Current  understanding  of  the  initiation  of 
detonation  of  heterogeneous  explosives  by  shock  therefore  divides  the  process  into  two 
distinct  stages: 

(i)  Ignition  of  a  small  fraction  of  the  explosive  at  random  sites  within  the  sample 

due  to  the  creation  of  hot  spots. 

(ii)  Growth  to  detonation  from  the  coalescence  of  the  energy  released  from  the 

individual  hot  spots. 

To  numerically  simulate  the  shock  initiation  of  a  condensed  phase  explosive  we 
therefore  need  to  find  appropriate  models  for  both  the  ignition  and  growth  stages,  and 
then  combine  these  into  an  appropriate  equation  for  the  overall  rate  of  explosive 
decomposition.  Many  models  for  hot  spot  formation  have  been  derived,  and  many 
different  reaction  rate  schemes  have  been  proposed.  A  critical  review  of  some  of  these 
reaction  rate  schemes  can  be  found  in  the  report  by  Jones  [45].  The  MULTI  code 
currently  contains  the  Forest  Fire,  Ignition  and  Growth,  and  CPEX  (Commercial 
Performance  of  Explosives)  reaction  rate  models,  although  additional  reaction  rate 
models  could  be  implemented  with  relative  ease. 

3.3.1  Forest  Fire 


The  Forest  Fire  model  was  one  of  the  first  reaction  rate  schemes  used  to  characterise 
the  decomposition  of  a  heterogeneous  explosive.  It  was  first  described  by  Mader  and 
Forest  [46],  and  a  good  description  of  the  model  can  be  found  in  the  monograph  by 
Mader  [5],  It  is  a  purely  phenomenological  model  and  assumes  that  the  reaction  rate 
can  be  expressed  as  a  power  series  in  pressure,  and  therefore  has  the  form: 


"  =  (i. 
dt 


f i=14 

^)exp  2  aiP1 

k  i=0  / 


(22) 


where  X  represents  the  fraction  of  reacted  explosive.  The  finite  rate  of  burning 
expressed  by  equation  (22)  yields  a  reaction  zone  of  finite  thickness  in  which  X  varies 
between  0  and  1,  and  which  contains  a  mixture  of  condensed  explosive  and  gaseous 
products.  In  such  regions  thermal  and  mechanical  equilibrium  between  the  two  phases 


6 


DSTO-TR-0705 


of  the  explosive  are  used  and  an  iterative  technique  is  used  to  determine  the  specific 
internal  energy  of  each  phase. 

Wedge  test  data  has  been  traditionally  used  to  obtain  the  coefficients  in  equation  (22). 
In  this  test  a  wedge  of  high  explosive  is  placed  in  front  of  a  planar  shock  wave 
generator  and  the  distance  that  the  shock  runs  through  the  wedge  before  detonation 
occurs  is  observed  with  a  streak  camera.  In  an  ideal  explosive  the  change  to  detonation 
is  marked  by  a  rapid  change  in  shock  velocity  and  is  easily  observed.  By  varying  the 
initial  shock  pressure  in  the  plane  wave  generator  a  plot  can  be  made  of  the  run 
distance  to  detonation  as  a  function  of  the  initial  sustained  shock  pressure.  Such  a  plot 
is  known  as  a  Pop  plot,  after  its  originator  A.  Popolato  [47].  The  Forest  Fire  coefficients 
are  then  obtained  by  fitting  to  the  Pop  plot  data.  In  order  to  perform  this  fit  several 
assumptions  have  to  be  made,  and  the  validity  of  these  assumptions  have  been 
analyzed  recently  by  Starkenberg  [48],  Lundstrum  [49],  and  Liang  et  al.  [50].  Although 
some  of  the  assumptions  regarding  the  flow  behind  the  shock  front  are  overly 
simplistic,  it  appears  that  replacing  these  assumptions  with  more  realistic  expressions 
has  little  effect  on  the  computed  coefficients.  Because  of  the  relative  ease  in  conducting 
the  wedge  test,  and  also  in  fitting  the  Forest  Fire  coefficients  to  the  Pop  plot  data,  it  is 
fairly  easy  to  obtain  the  Forest  Fire  coefficients  for  a  number  of  ideal  military 
explosives.  Mader  [5]  contains  a  comprehensive  listing  of  these. 

3.3.2  Ignition  and  Growth 

The  original  Lee  and  Tarver  model  [51]  divided  the  initiation  process  into  two  distinct 
stages.  In  the  first  stage,  the  ignition  stage,  the  passage  of  the  shock  front  creates 
hotspots  at  density  discontinuities  within  the  heterogeneous  material.  In  the  second 
stage,  the  growth  stage,  the  hot  spots  are  assumed  to  grow  by  a  grain  burning  process 
until  they  eventually  coalesce  to  form  a  stable  detonation.  The  model  is 
phenomenological  in  the  sense  that  plausible  assumptions  are  made  regarding  the 
physical  mechanisms  for  each  of  these  stages,  and  then  a  generalised  energy  release 
rate  of  the  following  form  is  considered: 

x ) 

-zf-  =  I(\-X)xr\r  +  G{\-X)X  pz ,  (23) 

ot 

r]=VJV,-\,  (24) 

where  X  is  the  fraction  of  explosive  that  has  reacted,  Vo  is  the  initial  specific  volume  of 
the  explosive,  Vj  is  the  specific  volume  of  the  shocked,  unreacted  explosive,  and  I,  x,  r, 
G,  y  and  z  are  constants. 

Different  hot  spot  models  lead  to  different  values  for  the  constant  r.  If  the  ignition  rate 
is  assumed  to  be  proportional  to  the  strain  rate  in  the  shocked  explosive  then  r  =  1.  If 
the  hot  spots  are  formed  by  the  stagnation  of  micro  jets  model  then  r  =  3,  while  if  the 
hot  spot  model  is  based  on  the  amount  of  plastic  work  done  during  void  collapse  then 
r  =  4.  Best  overall  agreement  with  experiment  has  been  found  using  r  =  4. 
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The  constants  x  and  y  are  related  to  the  choice  of  the  geometry  for  the  hot  spot 
combustion  process.  Hot  spots  can  be  considered  to  bum  outwards  from  the  void 
centre,  or  inward  over  the  total  grain  surface.  Lee  and  Tarver  considered  a  spherical 
hot  spot  burning  outward,  which  corresponds  to  y  =  2/3.  Requiring  that  the  rate  be  a 
maximum  when  the  combustion  surfaces  overlap  leads  to  the  value  x  =  2/9. 

The  remaining  constants  I,  G  and  z  are  found  by  combining  equation  (23)  with  the  JWL 
equation  of  state  (for  both  the  reacted  and  unreacted  explosive),  implementing  die 
scheme  in  a  Lagrangian  hydrocode  (DYNA2D),  and  then  fitting  simulated  pressure¬ 
time  records  to  experimentally  measured  embedded  gage  records.  In  their  original 
publication  Lee  and  Tarver  applied  the  model  to  the  explosives  PBX-9404,  TATB,  PETN 
and  cast  TNT,  and  found  that  they  were  able  to  obtain  good  agreement  with  a  variety 
of  experimental  data  for  explosives  subjected  to  sustained  shock  loading.  For  short 
pulse  duration  shock  initiation  experiments  however  the  model  needs  to  be  modified 
by  the  addition  of  an  extra  growth  term  before  good  agreement  with  experiment  is 
obtained.  Only  the  original  two  term  Ignition  and  Growth  model  is  implemented  in 
MULTI,  but  the  three  term  model  is  explained  in  Tarver  et  al.  [52],  and  could  easily  be 
added  to  the  code. 

Embedded  gauge  experiments  are  fairly  costly  to  conduct  however  and  so  relatively 
few  explosives  have  been  fitted  to  this  model.  Murphy  et  al.  [53]  have  addressed  this 
question  and  shown  that  it  is  possible  to  find  ignition  and  growth  parameters  for 
Composition  B  by  using  a  combination  of  existing  data  from  standard  tests  for  the 
material  and  extrapolation  of  the  remaining  unknown  parameters  from  similar  known 
explosives.  Simulations  of  the  wedge  test  and  failure  diameter  tests  were  found  to  be 
sufficient  to  define  the  ignition  and  growth  parameters  used  in  the  two  term  version  of 
the  reaction  rate  model.  The  coefficients  were  then  used  to  model  the  response  of 
several  two-dimensional  Composition  B  impact  initiation  experiments,  and  good 
agreement  with  the  experimental  data  was  found. 

Miller  [54]  has  also  considered  the  problem  of  determining  the  reactive  rate  parameters 
for  the  ignition  and  growth  model.  His  simplified  ignition  and  growth  model  (SIG) 
consists  of  only  two  adjustable  parameters,  the  ignition  (I)  and  growth  (G)  rate 
constants,  which  are  determined  from  experimental  data  on  failure  diameter  and  gap 
test  sensitivity.  Miller  has  applied  his  SIG  model  to  four  quite  different  explosives, 
PBXN-110,  PBXN-111,  PBXW-126  and  PBX  9404,  and  found  that  there  was  very  good 
agreement  between  simulated  and  observed  embedded  gauge  stress-time  profiles  for 
each  of  these  explosives. 

3.3.3  Programmed  Burn 

A  programmed  bum  is  used  in  a  hydrocode  to  simulate  the  detonation  of  an  explosive 
when  the  properties  of  that  particular  explosive  are  well  known,  and  it's  behaviour  is 
not  the  focus  of  the  investigation.  For  example,  a  programmed  bum  could  be  used  to 
detonate  a  charge  in  an  air  burst  calculation  when  the  effect  of  the  air  shock  on 
specified  targets  was  the  main  concern.  In  this  particular  application  the  characteristics 
of  the  explosive  are  well  known  and  the  programmed  bum  provides  a  simple  and 
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robust  method  of  ensuring  the  correct  detonation  pressure  at  the  explosive/air 
interface. 

In  a  programmed  bum  the  basic  assumption  is  that  the  detonation  wave  front  travels 
in  all  directions  at  the  Chapman-Jouguet  detonation  velocity.  Information  concerning 
the  energy  released  from  the  explosive  such  as  the  bum  time  (BT)  and  bum  interval 
(BI)  are  precalculated  and  stored  in  the  code  at  time  t  =  0.  For  example,  if  the  initiation 

point  of  the  detonation  is  at  xa,  then  the  arrival  time  of  the  detonation  at  a  typical 
computational  cell  centre  xi+i/2  is  given  by 

BT(i)  =  |xi+i  /2  -  Xd|  /  D  (25) 

where  D  is  the  detonation  velocity.  During  the  calculation,  when  the  cycle  time  tn  is 
greater  than  the  bum  time  BT  for  a  particular  cell,  but  less  than  the  time  (BT  +  BI),  then 
a  fraction  of  the  specific  energy  for  that  particular  explosive  is  deposited  into  the  cell. 
On  each  successive  cycle  this  process  is  repeated  until  the  cycle  time  is  greater  than  (BT 

+  BI).  The  fraction  of  specific  energy  deposited  each  time  is  given  by  (tn+1  -  tn)/BI,  and 
the  bum  interval  BI  is  chosen  so  that  the  reaction  zone  is  spread  over  several 
computational  cells.  In  this  way  the  correct  C-J  pressure,  density  and  velocity  can  be 
imposed,  and  the  resulting  detonation  front  is  guaranteed  to  have  the  correct  profile. 


4.  NON-IDEAL  EXPLOSIVES 

4.1  Introduction 

All  explosives  exhibit  some  form  of  diameter  effect  in  which  the  detonation  velocity 
decreases  with  decreasing  charge  diameter  until  a  critical  diameter  is  reached  and 
stable  detonation  cannot  be  sustained.  For  ideal  explosives  the  change  in  velocity  with 
diameter  is  minimal  until  very  close  to  the  critical  diameter,  whereas  for  non-ideal 
explosives  the  diameter  dependence  is  very  pronounced  and  the  velocity  at  the  critical 
diameter  can  be  as  low  as  30%  of  the  value  at  infinite  diameter.  Non-ideal  explosives 
contain  separate  fuel  and  oxidizer  species,  often  in  physically  separated  phases,  and 
their  heterogeneous  nature  leads  to  much  larger  reaction  zone  lengths  and  more 
curved  detonation  fronts  than  those  in  ideal  explosives. 

The  analysis  of  non-ideal  explosives  is  considerably  more  involved  than  that  of  ideal 
explosives,  and  has  only  begun  to  be  addressed  in  the  context  of  military  explosives 
during  the  last  few  years.  Reaction  rate  schemes  such  as  Forest  Fire  and  Ignition  and 
Growth  were  designed  for  ideal  military  explosives,  and  have  been  very  successful  in 
modelling  the  performance  of  these  types  of  explosives.  Their  applicability  to  non-ideal 
explosives  is  still  uncertain  however,  although  recent  progress  in  this  area  has  been 
made  by  Miller  [54],  and  Miller  and  Sutherland  [55]. 

Underwater  explosives  typically  display  highly  non-ideal  behaviour,  and  one  such 
explosive  which  has  been  of  considerable  interest  in  recent  years  is  PBXN-111 


9 


DSTO-TR-0705 


(formerly  known  as  PBXW-115),  which  has  the  formulation  20%  RDX,  43%  ammonium 
perchlorate,  25%  aluminium,  and  12%  HTPB  binder.  Kennedy  and  Jones  [56]  have 
analysed  the  detonics  of  this  underwater  explosive  using  small  divergent  detonation 
theory,  and  by  calibrating  a  reaction  rate  model  to  experimental  detonation  velocity 
measurements  as  a  function  of  charge  diameter.  The  model  was  implemented  in  the 
finite  element  hydrocode  DYNA2D,  which  was  then  used  to  simulate  a  variety  of 
initiation  and  detonation  tests,  and  the  results  were  generally  in  excellent  agreement 
with  the  experimental  data.  The  model  has  since  been  included  in  the  MULTI  code, 
and  the  next  section  briefly  describes  the  model  as  implemented  in  MULTI. 

4.2  Reaction  rate  scheme 

The  reaction  rate  used  to  model  PBXN-111  was  developed  for  composite  explosives  by 
Kirby  and  Leiper  [57]  and  has  the  form 


iX=a-x)\M  +  pa'  +  pa'l 

dt  L  T„  T;  Xf  _ 

(26) 

where 

(  3  y 

Phs  =  PA  P  for  p<  4pcr/3 

4  4d 

V  Fcr  ; 

(27) 

and 

phs  =  p  -  per  for  p>4pCr/3 

(28) 

X  is  the  progress  variable  which  describes  the  extent  of  reaction,  and  varies  from  0  for 
unreacted  explosive  to  1  for  detonation  products.  The  subscripts  are:  h  for  hotspot,  i  for 
intermediate,  and  f  for  the  final  stages  of  the  reaction.  There  are  four  adjustable 
parameters;  three  characteristic  reaction  times  t,  and  the  critical  pressure  pCr  that 
inhibits  the  onset  of  the  hotspot  reaction.  The  a;  are  Gaussian  form  functions  which 
describe  the  assumed  geometry  of  the  bum  front  and  which  switch  the  reactions  on 
and  off  as  the  various  phases  are  ignited  and  consumed.  The  a;  depend  on  X  and  the 
initial  formulation  of  the  explosive  as  follows: 

ah  =  exp{-[(X,-G)  / Wh]2}  for  0  <  X  <  Ci 

ah  =  exp{-[(^-Ci)/Wi]2}  for  Ci  <  X  <  1  (29) 

ai  =  0  for  0  <  X  <  G 

ai  =  1-  ah  for  Ci  <  X  <  Cf 

ai  =  exp[-[(?i-Ci)/Wh]2}  for  Cf  <  X  <  1  (30) 

af  =  0  for  0  <  X  <  Cf 

af  =  1-  ah  -ai  for  Cf  <  X  <  1  (31) 
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The  Gaussian  parameters  are  defined  in  terms  of  the  mass  fractions,  <E>,  of  the  three 
stages  as: 


G  =  <th2 

Cf  =  (Oh  +  <f>i)2  (32) 

Wh  =  2Q/(l  +  Ci) 

Wi  =  Oh(l-Oh) 

Wf  =  Of  (1  -  <J)f)  (33) 

The  equation  of  state  of  both  the  unreacted  and  reacted  phases  has  a  simple  density 
dependent  polytropic  equation  of  the  form 

e  =  -(g-\y  (34) 

P 

where  g  =  go  +  gap  +  g2P2  (35) 

The  values  of  the  constants  gi  for  the  explosive  products  are  found  by  fitting  to 
isentrope  data  from  an  ideal  thermodynamic  code,  while  the  constants  for  the 
unreacted  phase  are  obtained  by  considering  the  shock  data  of  the  explosive  and  using 
known  Hugoniot  data.  The  unreacted  solid  state  is  usually  a  mixture  of  ingredients, 
and  in  this  case  the  Hugoniots  of  the  components  are  combined  using  the  method  of 
Afanasenkov  et  al.  [58].  The  equation  of  state  for  the  partially  reacted  mixture  is  then 
completely  specified  by  invoking  the  mixture  rules 

P  =  Px  =  pP,  v  =  vx  =  vp,  e  =  (1-A,)ex  +  Xep  (36) 

where  the  subscript  x  refers  to  the  unreacted  explosive,  and  the  subscript  p  refers  to  the 
reacted,  or  product,  phase. 

5.  NUMERICAL  SOLUTION 


5.1  Introduction 

We  have  experimented  with  two  slightly  different  approaches  to  the  solution  of  the 
coupled  equations  described  in  the  previous  sections.  In  both  cases  operator  splitting 
was  used  to  decouple  the  transport  stage  from  the  chemical  reaction  stage,  and  then  the 
two-dimensional  transport  equations  were  further  decoupled  into  two  sets  of  one¬ 
dimensional  equations.  This  is  a  standard  technique  for  the  solution  of  coupled 
transport  equations  and  is  described  in  detail  by  Oran  and  Boris  [16]. 

In  the  first  version  of  the  code  the  density  of  each  material  on  the  grid  was  convected 
independently  by  making  repeated  calls  to  the  LCPFCT  algorithm  [59]  for  each 
material.  The  momentum  flux  and  total  energy  flux  (ie.  specific  internal  energy  plus 
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kinetic  energy)  for  each  material  were  then  added  and  the  combined  fluxes  convected 
using  LCPFCT.  This  approach  is  illustrated  in  the  equations  listed  below,  which  show 
the  sequence  of  equations  to  be  solved  for  the  x-pass  in  a  2D  cartesian  grid  calculation; 

Pi  is  the  density  of  the  i'th  material  on  the  grid,  pu  and  pv  are  the  x  and  y  components 
of  the  TOTAL  momentum  in  each  cell,  and  E  is  the  TOTAL  energy  per  unit  volume  in 
each  cell. 


5Pi  ,  5Piu  _  Q 
dt  dx 

Spu  |  0puu  0P 
dt  dx  dx 

(37) 

(38) 

dpv  |  0pvu  _  Q 

(39) 

dt  dx 

5E  5Eu  5Pu 

+  — 

dt  dx  dx 

(40) 

Because  of  the  multimaterial  capability  of  the  code  and  the  diffusive  nature  of  all 
numerical  transport  schemes  an  interface  tracking  algorithm  is  required  to  maintain  a 
sharp  interface  between  different  materials  on  the  grid.  In  the  first  version  of  the  code 
we  used  the  Simple  Line  Interface  Calculation  (SLIC)  scheme  of  Noh  and  Woodward 
[60].  SLIC  is  a  one-dimensional  alternating  direction  method  for  the  geometric 
approximation  of  fluid  interfaces.  The  scheme  constructs  a  representation  of  the 
interface  between  two  materials  from  the  volume  fractions  of  the  materials  in  the 
mixed  cell,  and  by  testing  whether  or  not  the  various  fluids  in  the  mixed  cell  are 
present  or  absent  in  the  zone  just  to  the  left  and  to  the  right  in  the  coordinate  direction 
under  consideration.  The  interface  between  two  different  materials  is  therefore 
represented  by  a  number  of  one-dimensional  components,  each  of  which  is  composed 
entirely  of  straight  lines  either  perpendicular  or  parallel  to  the  coordinate  direction 
under  consideration.  LCPFCT  assumes  that  the  different  materials  in  a  mixed  cell  are 
spread  homogeneously  throughout  the  cell,  and  hence  would  calculate  an  incorrect 
flux  in  a  mixed  cell.  The  SLIC  algorithm  determines  the  correct  amount  of  material  to 
be  advected  into  and  out  of  a  mixed  cell,  and  then  the  LCPFCT  algorithm  is  modified 
by  imposing  a  multiplicative  area  factor  at  the  correct  cell  boundary  for  each  material. 
Further  details  of  this  procedure  can  be  found  in  the  report  by  Milne  and  Carnegie  [32]. 

Whilst  the  above  scheme  convects  the  density  of  each  material  on  the  grid 
independently,  only  the  total  energy  in  each  cell  is  convected.  In  order  to  calculate  a 
unique  pressure  in  each  mixed  cell  both  the  density  and  specific  internal  energy  of  each 
material  in  the  cell  must  be  known,  and  so  a  way  of  dividing  the  total  internal  energy 
in  each  cell  between  the  component  materials  must  be  devised.  To  do  this  we  assumed 
pressure  equilibrium  between  the  materials  in  the  mixed  cell  (this  is  a  common 
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assumption  in  multimaterial  Eulerian  codes)  and  used  an  iterative  procedure  to 
determine  the  internal  energy  of  each  component. 

The  above  scheme  was  used  in  a  variety  of  trial  calculations,  including  both  one¬ 
dimensional  and  two-dimension  flying  plate  impact  simulations,  but  problems  were 
found  to  occur  at  the  interface  between  different  materials.  In  particular,  whilst  the 
average  numerical  values  for  the  dynamical  variables  were  accurately  predicted,  there 
were  perturbations  to  both  the  velocity  and  pressure  profiles  in  the  vicinity  of  the 
interfaces. 

These  disturbances  were  traced  to  the  way  in  which  the  total  energy  was  convected, 
and  to  a  discrepancy  between  the  mass  fluxes  and  interface  velocities  used  in  the 
partial  density  calculations,  and  the  momentum  fluxes  used  in  the  principle 
momentum  component  equation.  These  problems  were  eHminated  by  adding 
additional  coding  to  convect  the  internal  energy  of  each  species  separately,  and  by 
slightly  rearranging  the  order  of  the  transport  equations.  The  advantage  of  following 
the  internal  energy  of  each  material  is  that  an  iterative  calculation  is  no  longer  needed 
in  each  mixed  cell,  and  the  pressure  in  the  cell  is  determined  from  a  simple  volume 
average  of  the  partial  pressures  within  the  cell.  To  ensure  correct  volume  fluxes  it  was 
also  decided  to  base  the  interface  tracking  calculations  on  a  donor  cell  scheme  rather 
than  an  FCT  scheme.  Hence  in  the  later  version  of  the  code  the  fluxes  for  energy  and 
density  in  a  mixed  cell  are  overwritten  by  a  donor  cell  calculation. 

The  sequence  of  equations  to  be  solved  for  the  x-pass  in  a  2D  cartesian  grid  calculation 
in  this  version  of  the  code  now  have  the  form  shown  below,  where  pi  and  e,  are  the 
density  and  specific  internal  energy  respectively  of  the  i'th  material  on  the  grid,  and  u 
and  v  are  the  x  and  y  components  of  the  velocity  in  each  cell. 


dp i  ,  5pi u  _  0 
dt  dx 

(41) 

du  duu  du  VP 

— ■  H - =  u - 

dt  dx  dx  p 

(42) 

dv  dvu  du 

- 1-  — —  =  v  — 

(43) 

dt  dx  dx 

de.-  dei:u  _  du 

dt  dx  dx 

(44) 

The  current  version  of  MULTI  also  uses  an  interface  tracking  algorithm  based  on 
Youngs  method  [61].  This  is  a  two-dimensional  algorithm  which  represents  interfaces 
more  accurately,  as  each  portion  of  the  interface  in  a  cell  is  represented  by  a  straight 
line.  Simulation  results  on  test  problems  using  this  new  version  of  the  code  were  found 
to  be  both  more  accurate  and  more  robust,  and  hence  only  this  version  of  the  coding  is 
described  in  detail  in  the  following  section. 
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5.2  Program  MULTI 

Program.  MULTI  consists  of  a  number  of  separate  files,  each  of  which  will  be  described 
briefly  here.  The  main  driving  program  is  contained  in  multi.f .  This  reads  the  input  file 
det.inp,  calls  subroutine  INITAL,  which  initialises  all  the  variables,  then  calls 
subroutines  xtran.f  and  ytran.f,  which  transport  mass,  momentum  and  energy  in  the  X 
and  Y  directions  respectively,  and  then  calls  subroutine  burn.f,  which  models  the 
decomposition  of  any  explosive  material  on  the  grid,  multi.f  also  includes  function 
SPEED,  and  subroutines  MINTIME  and  TIMESTEP,  which  allow  the  user  to  specify  a 
variable  time  step  for  materials  described  by  the  HOM  equation  of  state. 

MULTI  allows  eight  different  material  types,  labelled  by  the  index  I,  and  characterised 
by  a  combination  of  several  different  equations  of  state  and  reaction  rate  laws.  The 
eight  different  material  types  are  listed  below: 


INDEX  I 
1 
2 

3 

4 

5 

6 

7 

8 


MATERIAL  TYPE _ 

Non  reactive  -  Perfect  Gas  equation  of  state 

Non  reactive  -  Mie-Gruneisen  equation  of  state 

Reactive  -  Forest  Fire  reaction  rate  model  and  HOM  equation  of  state 

Reactive  -  Ignition  and  Growth  reaction  rate  model  and  JWL  equation  of  state 

Reactive  -  an  early  version  of  JWL,  now  available  for  experimentation 

Reactive  -  Programmed  Burn  reaction  rate  model  and  JWL  equation  of  state 

Non  reactive  -  TAIT  equation  of  state  for  water 

Reactive  -  -  CPEX  burn  model  for  non  ideal  explosives 


multi.f  first  reads  the  input  data  file  det.inp  to  determine  the  number  of  materials  on 
the  grid  (NMAT),  the  type  of  boundary  conditions  (periodic  or  non-periodic),  and  the 
geometry,  which  is  specified  by  ALPHAX  and  ALPHAY.  ALPHA  =  1,  2  or  3 
corresponds  to  planar,  cylindrical,  or  spherical  geometry,  respectively,  and  ALPHAX 
and  ALPHAY  refer  to  the  geometry  in  the  X  and  Y  directions.  If  the  boundary 
conditions  are  non-periodic  then  they  can  be  either  reflective  (R)  or  transmissive  (T), 
and  the  boundary  conditions  at  the  extremes  of  the  X  and  Y  axes  are  specified  by 
setting  BCL,  BCR,  BCB  and  BCT  to  either  R  or  T.  The  time  step  can  be  either  fixed  or 
variable,  and  this  is  specified  by  setting  the  logical  variable  FIXED_TIME  to  be  either 
.true,  or  .false..  A  variable  time  step  is  currently  only  allowed  if  the  materials  on  the 
grid  are  described  by  either  the  perfect  gas  or  HOM  equations  of  state.  DELTAX  and 
DELTAY  specify  the  spatial  cell  sizes  in  the  X  and  Y  directions,  and  NCX  and  NCY 
specify  the  total  number  of  cells  in  the  X  and  Y  directions. 

det.inp  then  lists  information  about  each  of  the  materials  on  the  grid,  starting  with  the 
default  material,  which  is  automatically  given  the  highest  material  number  (equal  to 
NMAT).  The  default  material  (for  example,  air,  in  an  air  blast  calculation)  is  placed 
over  the  entire  grid,  and  then  other  materials  are  placed  at  specific  locations  on  the  grid 
by  specifying  their  minimum  and  maximum  X  and  Y  values  (IMIN,  IMAX  and  JMIN, 
JMAX).  Each  material  requires  specification  of  material  type,  initial  density,  initial 
velocity  components  Uo  and  Vo  (for  X  and  Y  directions  respectively)  and  PRIORITY. 
The  priority  value  determines  the  order  in  which  materials  are  convected  from  a  mixed 
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cell.  If  the  material  is  an  explosive  then  the  EMITTYPE  will  be  either  1  or  2,  depending 
on  whether  a  block  of  explosive  (rNTTTYPE=l)  or  circle  of  explosive  (INITTYPE=2)  is 
to  be  placed  on  the  grid.  If  a  block  of  explosive  is  used  then  NDUMPX  and  NDUMPY 
specify  a  subsection  of  the  explosive  in  which  an  additional  amount  of  internal  energy 
is  added  to  initiate  detonation.  If  a  circular  geometry  is  used  then  subroutine 
EXPLOSIVE_CIRCLE  (file  mkcirc.f)  works  out  (from  RADIUS,  XCENTRE,  and 
YCENTRE,  which  are  specified  in  det.inp)  which  cells  are  full  of  explosive,  and  which 
cells  have  fractional  amounts  of  explosive.  Subroutine  PUT_EXP_CIRCLE  (file 
put_circle.f)  then  fills  these  cells  with  appropriate  variable  values.  The  circular 
explosive  charge  is  initiated  by  compressing  a  smaller  circular  region.  The  radius  of 
this  region  is  typically  half  the  radius  of  the  full  charge,  and  the  explosive  is 
compressed  by  calling  subroutine  PUT_HP_CIRCLE  (file  put_hp_circle.f)  If  a 
programmed  burn  is  used  instead  to  initiate  the  explosive  then  the  compressive  radius 
is  set  to  zero. 

The  interface  tracking  part  of  the  code  is  contained  in  files  fracvol.f,  interfac.f,  and 
choke.f.  fracvol.f  is  the  basic  interface  tracking  routine  based  on  Youngs  method  [61]. 
interfac.f  is  used  by  fracvol.f  to  calculate  mesh  cell  wall /interface  intersections,  and 
choke.f  is  an  area  factor  which  controls  the  amount  of  material  advected  across  each 
cell  interface.  More  information  on  these  routines  can  be  found  in  the  report  by 
Carnegie  [33].  The  file  press.f  contains  the  equation  of  state  information  and  is  called  to 
calculate  the  pressure  for  each  of  the  material  types.  Some  of  the  equation  of  state 
coding  is  contained  in  files  jwl.f  and  homl.f,  which  are  self-explanatory,  homl.f 
contains  a  simple  addition  to  the  HOM  equation  of  state  to  ensure  consistency  with  the 
perfect  gas  equation  of  state  at  low  pressures,  burn.f  controls  explosive  burn.  It 
contains  subroutine  CP_BURN,  which  describes  the  reactive  rate  law  for  non  ideal 
explosives,  and  also  calls  ffire.f  and  igburn.f,  which  perform  Mader's  Forest  Fire  burn 
and  the  Lee  and  Tarver  Ignition  and  Growth  bum  respectively,  lcpdc.f  is  a  version  of 
lcpfct.f  which  solves  partial  density  and  energy  equations  assuming  a  donor  cell 
scheme  for  mixed  cells,  rtbis.f  is  a  real  function  which  calculates  roots  of  functions 
using  the  bisection  method,  and  is  used  in  several  subroutines  where  pressure 
equilibrium  is  assumed. 

5.3  Time  Step  Control 

MULTI  is  an  explicit  finite  difference  code  and  hence  the  time  step  8t  is  limited  by  the 
Courant  condition,  which  in  one  space  dimension  has  the  form  [16] 

8t  <  min  (  8x/{|v|  +  c})  (45) 

where  v  is  the  particle  velocity,  c  is  the  sound  speed,  and  8x  is  the  spatial  cell  size.  Use 
of  equation  (45)  ensures  the  stability  of  the  coupled  set  of  equations  (1)  -  (3),  (5)  and  (9). 
Expressions  for  the  sound  speed  can  be  derived  using  standard  thermodynamic 
arguments.  For  materials  described  by  the  FIOM  equation  of  state  Guirguis  and  Oran 
[62]  have  derived  expressions  for  the  sound  speed  in  the  condensed  phase,  gas  phase, 
and  a  mixture  of  both  phases,  and  these  expressions  have  been  coded  into  MULTI  so 
that  simulations  using  only  the  HOM  and  perfect  gas  equations  of  state  can  be  run 
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using  a  variable  time  step  calculated  from  equation  (45).  If  other  equations  of  state  are 
employed  in  the  calculation  then  a  fixed  time  step  must  be  used. 


6.  REPRESENTATIVE  RESULTS 

6.1  Bullet  impact  simulations 

MULTI  has  previously  been  used  to  simulate  the  impact  sensitivity  of  the  explosive 
Composition  B,  and  detailed  results  can  be  found  in  the  report  by  Borg  and  Jones  [35]. 
The  sensitivity  tests  used  steel  cylinders  of  varying  diameters  fired  against 
Composition  B  to  determine  the  critical  impact  velocity  for  initiation  as  a  function  of 
the  projectile  diameter.  The  experiments  were  later  modelled  by  Starkenberg  et  al.[63] 
using  Mader's  2DE  code,  and  good  agreement  was  found  between  the  calculations  and 
the  experimental  results.  The  MULTI  calculations  were  undertaken  partly  as  a 
validation  exercise  for  the  code.  As  the  simulations  were  performed  using  Mader's 
HOM  equation  of  state  and  the  Forest  Fire  reaction  rate  scheme  the  MULTI  results 
were  expected  to  agree  with  Starkenberg's  calculations,  apart  from  any  minor 
differences  arising  from  the  different  transport  algorithms  and  interface  tracking 
routines  in  the  codes. 

The  projectile  impact  calculations  were  performed  using  2D  cylindrical  geometry  and 
simulated  the  impact  of  a  the  cylindrical  steel  projectile  onto  a  cylindrical  Composition 
B  explosive  charge.  Both  the  projectile  and  explosive  charge  were  surrounded  by  air  at 
atmospheric  pressure.  The  air  was  described  by  a  perfect  gas  equation  of  state,  and  the 
steel  projectile  by  the  Mie-Gruneisen  equation  of  state. 

The  experiments  determined  the  critical  velocity  (Vc),  which  is  defined  as  the  velocity 
of  the  projectile  above  which  impact  results  in  shock  initiation,  and  below  which  the 
impact  results  in  a  lower  order  event  such  as  a  burn  or  deflagration.  The  critical 
velocity  is  a  function  of  projectile  diameter,  and  increases  as  the  diameter  of  the 
projectile  is  reduced.  Figure  1.  shows  a  comparison  between  the  simulated  results 
calculated  using  MULTI,  the  experimental  results  of  Slade  and  Dewey  [64],  and  the 
2DE  computations  of  Starkenberg  et  al.  [63].  The  MULTI  results  are  in  close  agreement 
with  the  calculations  of  Starkenberg  et  al.,  as  expected.  The  discrepancy  at  larger 
projectile  diameters  can  be  explained  by  differences  in  the  methods  used  to 
differentiate  between  a  detonation  and  a  fail  in  the  two  sets  of  calculations. 
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Figure  1:  Plot  of  critical  velocity  as  a  function  of  inverse  charge  diameter  for  the  MULTI 
simulations,  the  2DE  computations,  and  the  experimental  results. 


6.2  Underwater  Sensitivity  Test  Simulation 

MULTI  has  recently  been  used  to  simulate  the  Underwater  Sensitivity  Test  and 
Modified  Gap  Test  for  the  explosive  H-6  as  part  of  a  collaborative  effort  between  the 
TTCP  countries  to  validate  national  hydrocodes  for  the  prediction  of  underwater 
sympathetic  detonation.  The  US  is  the  lead  nation  in  this  activity  and  proposed  the 
numerical  simulation  of  the  underwater  sensitivity  of  H-6  as  a  common  benchmark 
problem,  and  has  supplied  each  nation  with  extensive  material  and  detonation 
property  data  on  H-6.  In  this  section  we  describe  some  of  the  recent  Underwater 
Sensitivity  Test  calculations  [37]  to  illustrate  the  application  of  MULTI  to  a  typical 
problem. 

The  Underwater  Sensitivity  Test  is  described  in  detail  by  Liddiard  and  Forbes  [65]. 
Figure  2  shows  a  schematic  of  the  test  assembly.  A  spherical  pentolite  donor  is 
suspended  in  the  aquarium  and  four  acceptors,  each  cylindrical  in  shape,  are  spaced  at 
varying  distances  around  the  donor  with  the  flat  surfaces  facing  the  donor  charge.  The 
donor  is  an  82.2  mm  diameter  sphere  of  cast  pentolite  (50%  TNT/50%  PETN)  weighing 
470  -  480  g.  The  acceptors  are  12.7  mm  thick  and  have  a  diameter  of  50.8  mm.  Chemical 
reaction  in  the  sample  explosive  is  detected  by  measuring  the  expansion  of  the  acceptor 
pellet  after  it  is  struck  by  the  shock.  If  the  axial  length  of  the  shocked  acceptor,  s,  is 
plotted  as  a  function  of  time,  t,  then  usually  a  fairly  constant  rate  of  expansion,  ds/ dt. 
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is  obtained.  The  slope  of  the  curve  is  then  used  to  characterise  the  degree  of  reaction  of 
the  sample  explosive. 

The  MULTI  code  was  used  to  simulate  the  UST  for  several  different  donor  -  acceptor 
gap  lengths.  The  simulations  were  rim  on  a  160  x  160  grid  in  2D  cylindrical  geometry 
with  8x  =  8y  =  0.1  cm  and  a  fixed  time  step  of  0.02  |ls.  The  detonation  of  the  pentolite 
donor  charge  was  modelled  using  a  programmed  burn  algorithm,  with  a  JWL  equation 
of  state]  to  describe  the  reaction  products.  For  pentolite  these  have  the  values  A  = 
5.4094,  B  =  0.093726,  Ri  =  4.5,  R2  =  1.1,  co  =  0.35,  and  Eo  =  0.081.  Initial  density  p0  =  1.70 
g/ cm3,  and  pq  =  0.255  Mbar.  The  pentolite  data  was  taken  from  Dobratz  and  Crawford 
[43].  The  water  was  modelled  by  the  Tait  equation  of  state,  as  described  in  section  2.2. 


TANK 

{60x60x60cm) 


472 -g 


Figure  2:  Schematic  illustration  of  the  Underwater  Sensitivity  Test. 


The  reactive  model  used  for  H-6  was  derived  by  McIntosh  and  has  been  described  in  a 
recent  publication  [66].  The  explosive  decomposition  is  modelled  by  the  Forest  Fire 
reaction  rate  law  and  was  calibrated  using  the  experimental  Pop  Plot  data  reported  by 
Hudson  [67],  McIntosh  used  the  JWL  equation  of  state  to  model  both  the  unreacted  and 
reacted  explosive  material,  and  the  DYNA  finite  element  code  was  used  to  model  both 
sensitivity  tests.  The  calculations  described  here  used  the  Forest  Fire  reaction  rate  law, 
but  the  HOM  equation  of  state  was  used  instead  of  JWL.  The  necessary  constants  were 
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provided  by  McIntosh  [66].  Table  1  shows  the  input  file  det.inp  for  a  typical  MULTI 
run.  In  this  case  the  donor  -  acceptor  separation  was  40  mm. 

Before  comparing  simulated  values  with  experimental  results  the  accuracy  of  MULTI 
with  respect  to  the  shock  pressure  in  water  was  checked.  Figure  3  shows  a  comparison 
between  the  simulated  pressure  and  the  experimental  pressure  from  the  Underwater 
Sensitivity  Test  calibration  reported  by  Liddiard  and  Forbes  [65].  There  is  quite  good 
agreement  between  the  MULTI  simulation  and  the  NSWC  calibration  for  all  the 
distances  of  interest. 


Pressure  Profile  in  Underwater  Sensitivity  Test 


Figure  3.  Shock  pressure  in  the  water  gap  of  the  UST  as  a  function  of  distance  from  the  surface 
of  the  pentolite  donor. 


Figures  4  and  5  show  plots  of  the  acceptor  thickness  as  a  function  of  time  for  water  gap 
distances  of  40  mm  and  55  mm  respectively.  These  show  a  fairly  constant  rate  of 
expansion,  and  this  is  also  found  experimentally.  The  simulated  ds/ dt  values  are  too 
high  though;  at  40  mm  and  55  mm  water  gap  distances  the  simulated  values  of  ds/dt 
are  1.39  mm/  ps  and  0.92  mm/  ps  respectively,  while  the  experimental  values  are  0.39 
mm/  ps  and  0.38  mm/  ps.  These  results  are  in  agreement  with  the  simulations  of 
McIntosh  however,  who  found  that  the  simulations  overestimated  the  velocities  by 
factors  of  about  two  to  three.  This  level  of  agreement  with  the  DYNA  results  is  to  be 
expected  as  basically  the  same  reactive  model  for  H-6  has  been  used  in  both  sets  of 
calculations. 
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The  rather  poor  level  of  agreement  between  the  simulated  results  and  the  experimental 
results  highlights  an  important  point  which  should  be  considered  in  any  reactive 
simulation;  namely,  the  level  of  agreement  will  only  be  as  good  as  the  suitability  of  the 
reactive  model  used.  Alternatively,  if  an  appropriate  reactive  model  is  used,  then  the 
model  should  be  calibrated  to  experimental  data  which  is  as  close  as  possible  in 
geometry,  shock  levels,  etc.,  to  the  experimental  data  which  is  being  simulated.  The 
Forest  Fire  model  for  H-6  was  calibrated  to  the  wedge  test  data  reported  by  Hudson 
[67].  Experimentally  the  wedge  test  has  almost  one-dimensional  flow,  while  the  UST  is 
very  much  a  two-dimensional  test.  Better  agreement  with  the  experimental  results 
would  have  been  obtained  if  the  reactive  model  had  been  calibrated  to  experimental 
data  exhibiting  a  strongly  two-dimensional  flow,  such  as  some  form  of  Gap  Test  data. 
This  point  is  discussed  in  more  detail  in  the  paper  by  Jones  [37]. 


Table  1:  Input  file  det.inpfor  a  typical  MULTI  run  for  the  Underwater  Sensitivity  Test. 

NMAT 

3 

PBC 

0.0 

ALPHAX 

1 

ALPHAY 

2 

BCL,  BCR,  BCB,  BCT 
■R’,  T,  ‘R’,  T 
DELTAX 
0.100 
DELTAY 
0.100 

DELTAT,  FIXEDJTIME 
0.0200,  .true. 

NCX 

160 

NCY 

160 

DEFEOS,  DEFC0,  DEFGAMMA,  DEFRHO,  DEFTMP 
7,  0.1483,  1.0,  1.0,  3.0E2 
DEFU0,  DEFV0,  DEFPRIORITY 
0.0,  0.0,  3 
INITTYPE,  RADIUS 
0,  0.0 

XCENTRE,  YCENTRE 

0.0,  0.0 

EOSTYPE(I),  C0(1),  GAMMA(I),  RHO0(1),  TMP0(1) 

5,  0.0,  0.0,  1.70,  3.0E2 
U0(1),  V0(1),  PRIORITY(I) 

0.0,  0.0,  1 

INITTYPE,  RADIUS 
2,  4.05 

XCENTRE,  YCENTRE 

0.0,  0.0 

IMIN(1),  IMAX(1),  JMIN(1),  JMAX(1) 
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0,  0,  0,  0 
ndumpx,  ndumpy 
0,  0 

EOSTYPE(I),  C0(1),  GAMMA(I),  RHOO(1),  TMP0(1) 
3,  0.0,  0.0,  1.760,  3.0E2 
U0(1),  V0(1),  PRIORITY(I) 

0.0,  0.0,  2 
INITTYPE,  RADIUS 
0,  0.0 

XCENTRE,  YCENTRE 

0.0,  0.0 

IMIN(1),  IMAX(1),  JMIN(1),  JMAX(1) 

81,  94,  1, 25 
ndumpx,  ndumpy 
0,  0 

MINSTEP 

1 

MAXSTEP 

3000 

IPRINT 

100 

PMIN 

0.0 

FULLDONOR 

.TRUE. 

DUMP,  RESTART 
.FALSE.,  .FALSE. 
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6.3  Modified  Gap  Test  Simulation 

This  section  briefly  describes  the  application  of  MULTI  to  the  simulation  of  the 
Modified  Gap  Test  for  the  TTCP  collaborative  hydrocode  validation  study  on  H-6 
described  in  the  previous  section.  The  Modified  Gap  Test  is  described  in  detail  by 
Liddiard  and  Forbes  [65],  and  Figure  6  shows  a  schematic  of  the  test  assembly.  It 
consists  of  a  cylindrical  donor  explosive,  a  PMMA  gap  of  variable  length,  and  an 
acceptor  explosive.  The  donor,  gap,  and  acceptor  are  all  50.8  mm  in  diameter,  and  the 
donor  for  the  H-6  acceptor  runs  was  50/50  pentolite.  Figure  6  shows  that  two  MGT  set¬ 
ups  are  fired  simultaneously  to  reduce  costs,  and  a  wooden  baffle  is  used  to  prevent 
gases  from  the  detonated  donors  from  obscuring  the  view.  The  effect  of  the  shock  from 
the  donor  explosive  on  the  acceptor  explosive  is  monitored  by  framing  camera 
observation  of  the  free  surface  behaviour  of  the  acceptor,  and  in  particular,  by 
measurement  of  the  axial  blow-off  velocity  of  the  acceptor. 

When  chemical  reaction  occurs  in  the  acceptor  the  free  surface  may  either  accelerate, 
decelerate,  or  stay  at  a  constant  velocity  over  the  entire  time,  depending  on  the  nature 
of  the  explosive  acceptor  and  the  degree  of  reaction.  The  free  surface  displacement  of 
the  acceptor,  xg,  is  measured  frame  by  frame,  and  then  the  free  surface  velocity  is 
obtained  by  measuring  the  slope,  dxg/dt,  over  a  linear  portion  of  the  xg-t  curve,  or  at  a 
point  on  the  curve  at  a  specified  value  of  xg  if  no  linear  portion  exists.  The  value  of  xg 
used  was  typically  somewhere  between  20  mm  and  40  mm,  but  unfortunately  does  not 
appear  to  be  specified  in  the  reports,  which  makes  detailed  comparison  with  the 
simulated  results  difficult. 

The  MGT  runs  for  H-6  were  modelled  using  a  Mie-Gruneisen  equation  of  state  for 
PMMA,  as  described  by  Bowman  et  al.  [68],  and  the  models  for  pentolite  and  H-6 
described  in  section  6.2.  The  simulations  were  run  in  2D  cylindrical  geometry  using  8x 
=  8y  =  0.1  cm  and  a  fixed  time  step  of  0.01  (is.  Grid  sizes  were  typically  of  the  order  of 
ncx  x  ncy  =  240  x  240.  The  H-6  acceptor  is  12.7  mm  thick,  the  pentolite  donor  is  65  mm 
thick,  and  for  H-6  the  experimental  results  used  PMMA  gap  lengths  varying  between 
32.81  mm  and  64.52  mm.  Table  2  shows  the  input  file  det.inp  for  a  MULTI  run  where 
the  PMMA  gap  had  a  length  of  35  mm. 

Figure  7  shows  the  Xg-t  curve  for  the  simulated  run,  and  Figure  8  shows  the  velocity¬ 
time  plot  obtained  by  numerically  differentiating  the  curve,  which  shows  that  Ua 
changes  continuously  with  time.  The  experimental  value  for  Ua  obtained  by 
interpolation  from  the  data  reported  by  Liddiard  and  Forbes  is  3.0  mm/gs,  but  the 
location  at  which  this  velocity  occurs  is  unknown.  McIntosh  [66]  has  modelled  the 
MGT  using  a  kinematic-isotropic  elastic-plastic  material  model  for  PMMA  and  used 
the  acceptor's  maximum  top  surface  velocity  for  comparison  with  the  experimental 
results.  If  we  adopt  this  approach  we  obtain  Ua  =  3.1  mm/gs  for  a  35  mm  gap,  which 
agrees  remarkably  well  with  the  experimental  value,  and  is  the  same  value  as  reported 
by  McIntosh.  For  longer  gap  lengths  the  agreement  is  less  encouraging,  and  reasons  for 
the  discrepancies  are  discussed  in  the  paper  by  Jones  [37]. 
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Table  2:  Input  file  det.inpfor  a  typical  MULTI  run  for  the  Modified  GapTest. 

NMAT 

4 

PBC 

0.0 

ALPHAX 

1 

ALPHAY 

2 

BCL,  BCR,  BCB,  BCT 
‘R’,  T,  ‘R\  T 
DELTAX 
0.100 
DELTAY 
0.100 

DELTAT,  FIXED.TIME 
0.0100,  .true. 

NCX 

160 

NCY 

160 

DEFEOS,  DEFCO,  DEFGAMMA,  DEFRHO,  DEFTMP 
1,0.03,  1.4,  1.0E-03,  3.0E2 
DEFUO,  DEFVO,  DEFPRIORITY 
0.0,  0.0,  4 
INITTYPE,  RADIUS 
0,  0.0 

XCENTRE,  YCENTRE 

0.0,  0.0 

EOSTYPE(I),  C0(1),  GAMMA(I),  RHOO(1),  TMP0(1) 

5,  0.0,  0.0,  1.70,  3.0E2 
U0(1),  V0(1),  PRIORITY(I) 

0.0,  0.0,  1 

INITTYPE,  RADIUS 
0,  0.00 

XCENTRE,  YCENTRE 

0.0,  0.0 

IMIN(1),  IMAX(1),  JMIN(1),  JMAX(1) 

1, 65,  1,  25 
ndumpx,  ndumpy 
0,0 

EOSTYPE(I),  C0(1),  GAMMA(I),  RHOO(1),  TMP0(1) 

2,  0.2432,  1.0,  1.180,  3.0E2 
U0(1),  V0(1),  PRIORITY(I) 

0.0,  0.0,  2 
INITTYPE,  RADIUS 
0,  0.0 

XCENTRE,  YCENTRE 

0.0,  0.0 

IMIN(1),  IMAX(1),  JMIN(1),  JMAX(1) 

66,  100,  1,  25 
ndumpx,  ndumpy 
0,0 
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EOSTYPE(I),  C0(1),  GAMMA(I),  RHOO(1),TMPO(1) 
3,  0.0,  0.0,1.760,  3.0E2 
U0(1),  V0(1),  PRIORITY(I) 

0.0,  0.0,  3 
INITTYPE,  RADIUS 
0,  0.0 

XCENTRE,  YCENTRE 

0.0,  0.0 

IMIN(1),  IMAX(1),  JMIN(1),  JMAX(1) 

101,  114,  1,  25 
ndumpx,  ndumpy 
0,  0 

MINSTEP 

1 

MAXSTEP 

6000 

IPRINT 

100 

PMIN 

0.0 

FULLDONOR 

.TRUE. 

DUMP,  RESTART 
.FALSE.,  .FALSE. 


Figure  6:  Schematic  illustration  of  the  Modified  Gap  Test  Assembly. 
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7.  DISCUSSION  AND  CONCLUSION 

During  the  time  that  MULTI  has  been  under  development  a  number  of  newer 
hydrocodes  have  become  available.  EFSAS  (Integrated  Fluid-Structure  Analysis  System) 
was  developed  by  Combustion  Dynamics  Ltd.  in  Canada  under  contract  to  Defence 
Research  Establishment  Suffield  [69].  EFSAS  began  as  a  pure  Eulerian  code  to  perform  air 
blast  calculations,  but  has  grown  considerably  since  then.  Currently  IFSAS  contains  a 
State  module  for  basic  gas  dynamics  calculations,  a  3D  CFD  module  for  compressible 
flows,  a  Blast  Raytracer  module  which  uses  image  charge  techniques  to  rapidly  compute 
pressure  loadings  for  internal  explosions,  and  a  finite  element  structural  response 
module  to  compute  elastic-plastic  deformation  of  stiffened  panels.  Additional  recent 
developments  include  implementation  of  a  fully  coupled  Euler-Lagrange  capability,  and 
multiphase  flow  algorithms  for  both  gaseous  and  condensed  phase  systems.  Current 
developments  include  implementation  of  detonation  and  shock  discontinuity  trackers, 
and  coupling  with  the  LSTC  DYNA3D  finite  element  code. 

CTH  is  an  Eulerian  code  developed  at  Sandia  National  Laboratories  to  model  three- 
dimensional,  multimaterial,  large  deformation,  strong  shock  wave  physics  [70].  CTH 
uses  tabular  or  analytic  equations  of  state  to  model  solid,  liquid,  vapour,  plasma,  and 
mixed-phase  materials.  Material  models  include  elastic-plastic  behaviour,  high 
explosives,  fracture,  and  motion  of  fragments  smaller  than  a  computational  cell.  The 
3D  version  of  CTH  uses  a  SLIC  interface  tracker,  while  the  2D  version  uses  a  high 
resolution  interface  tracker  which  assumes  that  the  interface  between  two  materials  can 
be  adequately  approximated  by  a  straight  line.  Explosives  are  described  by  a  JWL 
equation  of  state  and  a  programmed  bum  model. 

CAST  (Computational  Applied  Science  and  Technology)  is  a  suit  of  computer 
programs  which  has  been  developed  as  a  collaborative  venture  by  Fluid  Gravity 
Engineering  Ltd.  and  the  Defence  Evaluation  and  Research  Agency  of  the  UK  Ministry 
of  Defence  [71].  The  basic  modules  provide  a  fully  2D/3D  multi-material  hydrocode 
capability  including  material  strength  and  explosive  modelling  which  is  centred 
around  the  Euler  codes  EDEN  and  GRIM.  These  codes  have  the  ability  to  model  a  wide 
range  of  structural  dynamic  problems,  including  impact  studies,  response  to  air  blast, 
and  response  to  underwater  explosions.The  CAST  constitutive  models  take  into 
account  a  variety  of  real  material  effects,  incuding  work  hardening,  thermal  softening, 
and  strain  rate  dependency.  Other  modules,  such  as  a  fully  reactive  flow  capability,  can 
also  be  added. 

MESA  is  a  3D  Eulerian  code  which  treats  hydrodynamic  flow  and  the  dynamic 
deformation  of  solid  materials  [72].  It  was  written  at  Los  Alamos  National  Laboratory 
and  was  designed  specifically  for  simulations  of  armour  and  anti-armour  systems.  It 
incorporates  several  of  the  standard  strength  models  which  take  into  account  both  the 
elastic  and  plastic  regions  of  the  stress-strain  relationship  of  the  materials.  Explosives 
are  described  using  a  programmed  bum  model.  MESA  has  been  used  to  study  the 
formation  of  non  axisymmetric  shaped  charge  jets,  the  penetration  of  reactive  armour 
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by  jets,  and  long  rod  penetrators.  MESA  contains  highly  accurate  interface  trackers 
which  can  resolve  highly  contorted  deformations  when  only  four  computational  cells 
are  used  to  define  plate  thickness. 

AUTODYN-2D  and  AUTODYN-3D  were  developed  by  Century  Dynamics  Ltd.  and 
are  interactive,  integrated  hydrocodes  which  provide  a  number  of  fully  coupled 
numerical  processors  [73].  These  include  Lagrangian,  Eulerian,  Smooth  Particle 
Hydrodynamics,  and  Arbitrary  Lagrange  Euler  processors.  The  codes  are  particularly 
suited  to  the  modelling  of  impact,  penetration,  blast  and  explosive  events.  The  Euler 
codes  use  a  version  of  the  SLIC  interface  tracker  to  define  material  interfaces,  and  an 
artificial  viscosity  to  spread  shock  fronts  over  several  mesh  intervals.  Many  material 
equations  of  state  and  material  strength  models  are  available  in  the  codes.  The  strength 
models  include  Hydrodynamic,  Elastic,  Piecewise  hardening,  Brittle,  and  also  the 
Johnson-Cook,  Zerilli-Armstrong,  Steinberg-Guinan  and  Johnson-Holmquist  models, 
which  include  the  effects  of  strain  and  strain  hardening  and  thermal  softening. 

WSD  staff  are  currently  using  IFSAS,  LS-DYNA3D,  and  MULTI  for  a  variety  of 
simulation  tasks  in  the  areas  of  warhead  simulation,  blast/target  interactions,  and 
explosives  simulation.  IFSAS  has  a  user  friendly  interface  and  is  ideal  to  use  when 
standard  shock/material  interaction  problems  need  to  be  modelled.  MULTI,  on  the 
other  hand,  was  designed  to  provide  an  in-house  code  which  could  be  quickly  adapted 
to  model  more  novel  explosive /shock  problems.  Recent  examples  within  WSD  have 
included  the  design  and  implementation  of  an  energetic  model  for  the  explosive 
PBXN-111  (formerly  PBXW-115)  [56],  and  fully  three-dimensional  air  blast  simulations 
[28].  MULTI  is  written  in  FORTRAN  and  runs  from  a  standard  input  data  file,  as  the 
examples  in  the  previous  section  have  shown,  but  users  are  expected  to  be  familiar 
with  CFD  techniques  and  to  be  capable  of  modifying  the  coding  when  non-standard 
problems  arise.  The  ability  to  modify  the  source  code  in  this  way  is  one  of  the 
advantages  of  MULTI  which  is  not  shared  by  the  more  commercial  codes  listed  above, 
although  some  of  them,  for  example  AUTODYN  and  CAST,  do  allow  users  to  write 
their  own  subroutines  to  specify  particular  material  models.  A  further  advantage  of  an 
in-house  code  such  as  MULTI  is  the  ability  to  tailor  the  code  to  suit  the  architecture  of 
whichever  computer  system  is  available  at  the  time.  For  example,  a  simplified  version 
of  MULTI  has  already  been  configured  to  run  in  parallel  mode  on  the  AMRL  POWER 
CHALLENGE  XV  [74].  Further  development  of  MULTI  will  depend  on  the  interests 
and  expertise  of  WSD  staff,  and  the  future  direction  of  WSD  tasks. 
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